function [llsT1Map, llsEquiliMagnMap, xVFAFitTerm, yVFAFitTerm] = fitLLSVFAdata(rawVFAdata, flipAngleMaps, acqVFAParam, TR)
%FITLLSVFA Calculates T1 from Monte Carlo VFA data from the linear 
%y = mx+b form of equation.
%   This function output the T1 values extracted from fitting the 

VFAdata = double(rawVFAdata);

xVFAFitTerm = VFAdata./tand(flipAngleMaps); % x = Signal/tan(alpha)
yVFAFitTerm = VFAdata./sind(flipAngleMaps); % y = Signal/sin(alpha)

% Fit for y = mx + b
% fitParamVFA(:,1) = m
% fitParamVFA(:,2) = b


fitParamVFA = zeros(acqVFAParam.xDim, acqVFAParam.yDim, acqVFAParam.zDim, 2);

for x = 1:acqVFAParam.xDim
    disp(x)
    for y = 1:acqVFAParam.yDim
        for z = 1:acqVFAParam.zDim
            fitParamVFA(x,y,z,1:2) = polyfit(squeeze(xVFAFitTerm(x,y,z,:)), squeeze(yVFAFitTerm(x,y,z,:)), 1);
        end
    end
end

longTR = zeros(acqVFAParam.xDim, acqVFAParam.yDim, acqVFAParam.zDim);
longTR = TR;
% m = exp(-TR/T1), solve for T1
llsT1Map = -longTR./log(fitParamVFA(:,:,:,1));

% b = M0(1-exp(-TR/T1)), solve for M0
llsEquiliMagnMap = fitParamVFA(:,:,:,2)./(1-exp(longTR./llsT1Map));

end

